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Abstract 

A new method for solving numerically stochastic partial differential equa- 
tions (SPDEs) with multiple scales is presented. The method combines 
a spectral method with the heterogeneous multiscale method (HMM) pre- 
sented in [W. E, D. Liu, and E. Vanden-Eijnden, Comm. Pure Appl. Math., 
58(11) :1544-1585, 2005]. The class of problems that we consider are SPDEs 
with quadratic nonlinearities that were studied in [D. Blomker, M. Hairer, 
and G. A. Pavliotis, Nonlinearity, 20(7):1721-1744, 2007.] For such SPDEs 
an amplitude equation which describes the effective dynamics at long time 
scales can be rigorously derived for both advective and diffusive time scales. 
Our method, based on micro and macro solvers, allows to capture numeri- 
cally the amplitude equation accurately at a cost independent of the small 
scales in the problem. Numerical experiments illustrate the behavior of the 
proposed method. 

Keywords: Stochastic Partial Differential Equations; Multiscale Methods; 
Averaging; Homogenization; Heterogeneous Multiscale Method (HMM) 



1. Introduction 



Many interesting phenomena in the physical sciences and in applications 
are characterized by their high dimensionality and the presence of many 
different spatial and temporal scales. Standard examples include atmosphere 
and ocean sciences 23|, molecular dynamics 16[ and materials science 17 
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The mathematical description of phenomena of this type quite often leads 
to infinite dimensional multiscale systems that are described by nonlinear 
evolution partial differential equations (PDEs) with multiple scales. 

Often physical systems are also subject to noise. This noise might be 
either due to thermal fluctuations 1^| , noise in some control parameter 18 



coarse-graining of a high-dimensional deterministic system with random ini- 



tial conditions [26|, |35| , or the stochastic parameterization of small scales [14 
High dimensional multiscale dynamical systems that are subject to noise can 
be modeled accurately using stochastic partial differential equations (SPDEs) 
with a multiscale structure. There are very few instances where SPDEs with 
multiple scales can be treated analytically. The goal of this paper is to 
develop numerical methods for solving accurately and efficiently multiscale 
SPDEs. Several numerical methods for SPDEs have been developed and an- 



alyzed in recent years, e.g. (4J, [13|, |31j, based on a finite difference scheme in 
both space and time. It is well known that explicit time discretization via 
standard methods (e.g., as the Euler-Maruyama method) leads to a time-step 
restriction due to the stiffness originating from the discretisation of the dif- 
fusion operator (e.g. the CFL condition At < C(Ax) 2 , where At and Ax are 
the time and space discretization, respectively). The situation is even worse 
for SPDEs with multiple scales (e.g. of the form (J2]) and fl3]) below) as in this 
case the Courant-Friedrichs-Lewy (CFL) condition becomes At < C(Ax-e) 2 , 
where e 1 is the parameter measuring scale separation. Standard explicit 
methods become useless for SPDEs with multiple scales. 

Such time-step restriction can in theory be removed by using implicit 



methods as was shown in [31|. However the implicitness of the numerical 
scheme forces one to solve potentially large linear algebraic problems at each 
time step. Furthermore, it was shown in 2l| that implicit methods are not 
suited for studying the long time dynamics of fast-slow stochastic systems as 
they do not capture the correct invariant measure of the system. Although 
this result has been obtained for finite dimensional stochastic systems, it 
is expected that it also applies to infinite dimensional fast-slow systems of 
stochastic differential equations (SDEs), rendering the applicability of im- 
plicit methods to SPDEs with multiple scales questionable. We also note 
that a new class of explicit methods, the S-ROCK methods, with much 
better stability properties than the Euler-Maruyama method was recently 
introduced in |l|, |2|, y|. Although these methods are much more efficient than 
traditional explicit methods, computing time issues will occur when trying 
to solve SPDEs with multiple scales as considered here, since the stiffness 
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is extremely severe for small e. Furthermore, capturing the correct invariant 
measure of the SPDE for At > e is still an issue for such solvers. 
In this paper we consider SPDEs of the form 



dtv = Av + F(v) + eQ£, 



(1) 



posed in a bounded domain of M. with appropriate boundary conditions. The 
differential operator A is assumed to be a non-positive self-adjoint operator 
in a Hilbert space H, £ denotes space-time Gaussian white noise, Q is the 
covariance operator of the noise and we take e C 1. We assume that the op- 
erator A has a finite dimensional kernel, M . This assumption leads to scale 
separation between the slow dynamics in M and the fast dynamics in the 
orthogonal complement of the null space Af ± , where % — J\f © M^. In this 
paper we will furthermore assume that noise acts directly only on the orthog- 
onal complement J\f . When noise acts also on A/", different distinguished 
limits than the ones considered in this paper should be considered. 

In order to describe the longtime behavior of the SPDEs We perform an 
advective rescaling set v(t) := eu(et). Using the scaling properties of white 
noise we obtain the following singularly perturbed SPDE 



Another scaling is of interest, namely the diffusive rescaling v(t) := eu(e 2 t) 
which leads to the SPDE 



For concreteness, we will focus on the class of SPDEs with quadratic nonlin- 
earities that was considered in 0], and assume that F{u) — f(u) + e a g(u), 
where / is a quadratic function (e.g. f(u) = B(u,u), a symmetric bilinear 
form), g a linear function and the exponent a is either 1 or 20 The choice of 



Usually the functions / and g involve derivatives of the function u. For example, 
for both the Burgers and the Kuramoto-Shivashinsky equation we have f(u) — ud x u. 
The linear function g(u) is included to induce a linear instability to the dynamics. In 
the case of the Burgers equation we will simply take g(u) = u whereas in the case of 
the Kuramoto-Shivashinsky equation we can take g(u) = d^u. Further discussion can be 
found in Section H and in [301 ] . 






(3) 
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a will depend on the particular scaling. Singularly perturbed SPDEs with 
quadratic nonlinearities provide a natural testbed for testing the applicabil- 
ity of the heterogeneous multiscale method to infinite dimensional stochas- 
tic systems, since a rigorous homogenization theory exists for this class of 
SPDEs [zj]. Furthermore, SPDEs of this form arise naturally in stochastic 



models for climate |23[ and in surface growth |19|, [34J . Finally, it has already 
been shown through rigorous analysis and numerical experiments that these 
systems exhibit a very rich dynamical behavior, such as noise-induced tran- 



sitions [33] and the possibility of stabilization of linearly unstable modes due 
to the interaction between the additive noise and the scale separation jsl, |io[ • 
We believe, however, that the methodology developed in this paper has a 
wider range of applicability and is not restricted to SPDEs with quadratic 
nonlinearities. Further comments about the class of SPDEs for which we 
believe that the proposed numerical method can be applied can be found in 
Section [5j 

Our numerical algorithm is based on a combination of a spectral method 
with micro-macro time integration schemes. We denote by x = P c u the 
projection onto M and by y = P s u, P s = I — P c the projection onto J\f . We 
then rewrite ([2]) and fl3]) as fast-slow system of SDEs 

x = a(x,y), (4a) 

y = -Ay + b(x,y) + - r QC, (4b) 
e v e 



and 



x = -a(x,y), (5a) 
y = \Ay + -b{x,y) + -Q£, (5b) 



e 2 e e 



where the functions a(x, y) and b(x, y) are the projections of F(u) onto M 
and A/" -1 . We remark that an 0(1) nonlinear term can be added in f[5al) . 
The fast-slow systems @j and dSJ) resemble fast slow systems for SDEs [29|, 
Ch. 10,11]. However, the fast process y is infinite dimensional and the well 
known averaging and homogenization theorems 0, [28[ do not apply. 

Averaging and homogenization results for SPDEs have been obtained 
recently |l(J, |7|. In particular, provided that the fast process y in (@J has 
suitable ergodic properties, then the slow variable x converges, in the limit 



4 



as e tends to 0, to the solution of the averaged equation 



x = a(x), (6) 

where the averaged coefficient is given by the average of a(x, y) with respect to 
the invariant measure of the (infinite dimensional) fast process y. When this 
average vanishes (i.e. the centering condition from homogenization theory is 
satisfied) then the dynamics at the advective time scale becomes trivial and it 
is necessary to look at the dynamics at the diffusive time scale, equations (JS}. 
It was shown in [tJ that the slow variable x of this system of equations, the 
solution of (15 ap , converges in the limit as e tends to to the solution of the 
homogenized equation 

± = a(x) + a(x)W, (7) 

with explicit formulas for the homogenized coefficients-see Section [3] for de- 
tails. For finite dimensional fast systems, the coefficients in fl^]) and (j7j) can be 



calculated, in principle, in terms of appropriate long-time averages-see [29 
for details. The numerical method proposed in 32] and analyzed in 
coined under the name of the Heterogeneous Multiscale Method (HMM), 
relies on the numerical approximation of the coefficients in ([6]) and flTj) by 
solving the original fine scale problem on time intervals of an intermediate 
time scale and use that data to evolve the slow variables using either 
or (171). In this paper we show how this methodology, when combined with a 
spectral method, can also be applied to SPDEs with multiple scales, that is, 
to the systems (j3J) and (jSJ). The aim of the present paper is to present the 
algorithm and report numerical experiments. The analysis of the proposed 
numerical method and the extension to more general classes of SPDEs with 
multiscale structure will be presented in a forthcoming paper. 

The rest of the paper is organized as follows. In Section [2] we present 
our new algorithm. Analytical and computational techniques for the analysis 
of SPDEs with multiple scales at the heart of the multiscale algorithm are 
presented in [3J In Section H] we present numerical experiments. Section is 
reserved for conclusions and discussion on future work. 



2. Numerical method 

We propose a numerical algorithm to approximate numerically the so- 
lution of (JT]) based on a micro-macro algorithm, capable of capturing the 
effective behavior of the SPDE. We explain the numerical algorithm for the 
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case of diffusive time scale (the hardest numerically) and comment on the 
advective time scale later in this section. 

2.1. Multiscale Algorithm 

We consider SPDEs ([I]) in a Hilbert space Ti with norm || ■ || and inner 
product (-,-). A denotes a differential operator, £ space-time white noise and 
Q the covariance operator of the noise. We assume that A is a self-adjoint 
nonpositive operator on % with compact resolvent. We denote its eigenvalues 
and (normalized) eigenfunctions by {—\k, ek}^ =1 : 

- Ae k = X k e k , k = 1, . . . (8) 

The eigenfunctions of A form an orthonormal basis in H. We assume that 
A and the covariance operator of the noise Q commute. Thus, we can write, 
formally, 

+oo 

fc=l 

where {£,k(t)}^=i are independent one-dimensional white noise processes, i.e., 

mean-zero Gaussian processes with ((,k(t)£,j(s)) = SkjS{t — s), k, j — 1, 2, 

Here 5k j and 5(t — s) are the usual Kronecker delta functions. 

Furthermore, we will assume that A has a finite dimensional kernel M := 
{h e% : Ah = 0}, dim(AT) = N < +oo and write U = M © Af ± . We 
introduce the projection operators 

V c : U^N (10a) 
V S = I-V C : H^N X : (10b) 

and write x := P c u, y := P s u. Finally, we will assume that noise acts only 

on A/^, i.e. q k = 0, k = 1 . . . N. 

Step 1. Decomposition in a fast-slow system. 

Using the projection operators defined in f llOap and (UObj) . equation ([T]) can 
be written as a fast-slow stochastic system 

x = -V c F(u), (11a) 
y = -^Ay + -V c F(u) + -QZ, (lib) 
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where x(t) e M. N since dim (A/") = N. We order the pairs of eigenfuctions and 
eigenvalues such that the kernel M is spanned by the first N eigenfunctions 
of A. We can write 

N +00 

x = ^x k e k and y = ^ Vk&k- 

k=l k=N+l 

Moreover, we introduce 

a k (x,y) := (V c F,e k ) for 1 < k < N, (12) 
b k (x jy ) := (V s F,e k ) iork>N. (13) 

Remark 2.1. As mentioned in the introduction we will often consider the 
case F(u) = f(u) + e 2 g(u). Then the above decomposition reads 

a k (x,y) := (PJ,e k ) + e{P c g,e k )=a k l {x,y) + ea k 1 {x,y), (14) 
b k (x,y) := (Vsf,e k ) + 6{V s g,e k )=b k (x,y) + eb k {x,y). (15) 

where we notice that for a linear function g(u) = vu we simply have a k (x, y) = 
vx k ,b k (x,y) = vy k . 

Then, in view of (jSJ) and we can rewrite the system (|lip in the form 

x k = -a k (x,y), k = l,...N, (16a) 
e 

11 1 

Vk = —^^kyk + -b k (x,y) + -q k ^ k , k = N + 1, N + 2, . . . (16b) 

Equations pip , resp. f fl6|) . are the infinite system of singularly perturbed 
SDEs that we want to solve numerically. 

Step 2. Truncation. 

We consider a finite dimensional truncation of the above system and keep M 
fast processes H : 

x = -a(x,y), (17a) 
e 

y = -vW+-b(x,y) + -Q M £, (17b) 



2 To simplify the notations we will use a new labeling of the index for the truncated fast 
system and write (y±, . . . , yA/)instead of (un+i, ■ • ■ , Dn+m) and similarly for the eigenval- 
ues Afe and the noise intensity qk- 
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where x = (x l7 . . . , x N ) T , y = (y u y M ) T , $, = . . . , £ M ) T and 

a(x,y) = (a 1 (x,y),...,a 7V (x,y)) T , (18) 
b(x,y) = (6 1 (x,y),...,6 A/ (x,y)) T , (19) 

and Am = diag(Ai, . . . , Am) and Q M = diag(gi, . . . , q^). For the decompo- 
sition f lT5|) . (TH|) . we will use the notations 



a(x,y) = a (x,y) + eax(x,y), (20) 
b(x,y) = b (x,y) + eb 1 (x,y), (21) 

where ao, ai G ~R N and bo, bi G M A/ with components similar as in ( TT8|) or 

Step 3. Numerical solution of the reduced system. 

The reduced system (fTTl) is solved by a micro-macro algorithm following 



32l . Il4j . It consists of a macrosolver (we use the notation X n := X(t n )), 



chosen here to be the Euler-Maruyama scheme 

X n+1 = X n + Ata M + <r n M AW n , (22) 

where AW n (the Wiener increment) is jV(0, At). Notice that At represents 
here a macrotime step, i.e., At can be chosen much larger than e. The drift 
function a M ~ a M (X n ) and diffusion function cr M ~ <r M (X n ) appearing in 
( 122]) . recovered from a time-ensemble average, are given by 



j k e T +L-i 



I fa K n-T+L-1 V 
+ ^7^~~ <9x a (^n, ^ 1 ,^',j) a (^ra) ^n/j)> 

(23a) 

if £ T +L-1 L' 



(23b) 



S 



where Y l ,Y 2 are the solutions of a suitable auxiliary system (given in 
below) involving the fast problem (117bj) . Here K denotes the number of sam- 
ples taken for the numerical calculation, L, V the number of micro timesteps 
and It a number of initial micro timesteps that are omitted in the averaging 
processes to reduce the effect of transients (see below). 

Auxiliary system. As observed in 32], for diffusive timescales, computing 
effective coefficients via time-averaging (relying on ergodicity), may require 
to solve ( 117bj) over time T = 0(e~ 2 ). To overcome this problem, it was 
suggested again in |32| to replace the fast process in ( 117bj) by (y ~ y 1 + ey 2 ) 



y 1 = -\^Aiy l + -Qm£, (24a) 



e 2 - e 



y 2 = -^AA/y' + ^b^y 1 ). (24b) 



The numerical approximations Y 1 ,Y 2 of (I24al) and f!24bl) . respectively, are 
the functions appearing in the averaging procedure to recover the macro- 
scopic drift and diffusion functions (see (I23ap )-( l23bl) ). Notice that we fix the 
slow variables in the system ( I24bl) at the current macro state X n . We use 
again the Euler-Maruyama method and compute Y 1 , Y 2 as 



YnMi = Y} e - S 4A M Y} e +^±Q M J n , (25a) 



e 2 ' e 



Y n,e+i - ^ _ ^"Aa./F^ + — b(X n , F n ^), (25b) 

where J„ = diag(J^, . . . , J* f ) and is a W(0, 1) random variable. The index 
n refers to the macrotime, t n . We compute f)25al) over L + L' microtime steps, 
(125bl) over L microtime steps to compute the time-ensemble average (I23ap . 
Notice that for the microsolver, the timestep 6t resolves the fine scale e 2 . The 
initial values for the micro solver are taken to be (for n > 1 

v 1 — v 1 v 2 — v 2 

I n,0 ~ 1 n-l,e T +L+L'-n 1 n,0 — 2 n-l,£ T +L-H 

and Yq Q = Yq = for n = 0. The motivation for computing the above time 
averages is given in the next section. 

Remark 2.2. We notice that the auxiliary system (12^1) is degenerate, since 
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the noise in (11 7b[) is additiveu This implies that the results presented in 14- 
App. B] are not applicable in this case and a more elaborate analysis is 
required for proving geometric ergodicity. This analysis, based on the ergodic 



theory for hypoelliptic diffusions }2a] , will be presented elsewhere. In the 



present paper we will assume that the auxiliary process 023]) is ergodic. 

Advective time scale.. A similar algorithm can be derived for the advective 
time scale. We consider the fast-slow system (jlj) that after projection and 
truncation reads 

x = a(x,y), (26a) 

y = -iA M y + b(x,y) + -Q M £, (26b) 
e e 

similarly as ( TT71) . The macrosolver, chosen to be the Euler explicit method, 
is given by 

X n+ i = X n + Ata^j, 
where the effective force is given by the time average 



al/ = E a(X n ,F n/j ), (27) 

3 = 1 t=tr 

(28) 

where Y n ij is a numerical approximation of the truncated fast system (135bl) 
with a slow variable fixed at time t n . As previously, K denotes the number of 
samples and L the number of micro timesteps and £t is the number of initial 
micro timestep ommited to reduce the transient effects. For the advective 



scaling, there is no need for an auxiliary problem for the micro solver [32 



3. Averaging and Homogenization for SPDEs 

In this section we summarize recent results on the averaging and homog- 
enization for SPDEs 10, 0] that are the analytical foundation on which the 



numerical algorithm presented in Section [2] is built. 



3 Indeed, the auxiliary system in [32|, [l4[ will always be degenerate, whenever the noise 



in the fast/slow system of SDEs that we want to solve is additive. 
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3.1. Analytic form of the homogenized coefficients 

In this section we briefly discuss the analytical form of the effective system 
corresponding to ffTTl) . Under the assumption that the vector field a (x, y) 
(see (1201) ) is centered with respect to the invariant measure of the fast process, 

/ a (x,y)/i(dy) = 0, (29) 
Jrm 

then the slow process converges to a homogenized equation of the form 

dX = a M (X) dt + & M (X) dW, (30) 

where W represent an iV— dimensional Wiener process and the SDE (1301) is 
interpreted in the Ito sense. The subscript M are used to emphasise the fact 
that the homogenized coefficients depend on the number of fast processes 
that we take into account. An analytic expression for the coefficients that 
appear in f[3"Ul) is given by 

a M (x) = lim / ^(rfy 1 ,rfy 2 )V,a(x,y 1 )y 2 (31a) 
lim / nidy 1 ) / E y i V x a(x, y\ 2 Ja(x, y 1 ) ds, 

e^O J RM J 



+ 



o- m (x)(<7m(x)) t = 2 lim / /i(rfy 1 )a(x,y 1 ) 

e^O J RM 
r+oo 

(8) / E y ia(x,y e 1 2 Jrfs. (31b) 
Jo 

Here fi(dyi) denotes the invariant measure of the process y 1 which is given 
by (I3"3"|) and ^(dy 1 , dy 2 ) denotes the invariant measure of the the process 
{y\ y 2 }- Notice that = y* is the solution of the rescaled process corre- 
sponding to f)24al) . i.e., y 1 = — A^y 1 + Qm£,- Alternatively, the calculation of 
the coefficients slm{x) and <tm{x) which appear in the homogenized equation 
can be obtained by the solution of the Poisson equation 

- C M (p = a (x,y), (32) 

where Cm is the generator of the fast (truncated) Ornstein-Uhlenbeck pro- 
cess. This process is ergodic and its invariant measure is Gaussian: 

K d y) = J_1 q2j d y> ( 33 ) 
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where Zj^ denotes the normalization constant. We notice that the sys- 
tem (TTTj) is a finite dimensional fast-slow system of SDEs for which stan- 



dard homogenization theory applies [5|, [28|, [29[ . For quadratic nonlinearities 
the Poisson equation (132]) can be solved analytically. The calculation of the 
coefficients in the homogenized (amplitude) equation reduces then to the cal- 
culation of Gaussian integrals that can also be done analytically. This will 
be done in Section I3~3"l 

3.2. The Advective Time Scale 

Averaging problems for fast-slow systems of SPDEs were studied recently 



in 



lOf] and their results can be applied to (jl]). One important observation is 
that in the system (jl]), the fast process is, to leading order C(l/e), an infinite 
dimensional Ornstein-Uhlenbeck process. The ergodic properties of such an 
infinite dimensional process can be analyzed in a quite straightforward way 
and the invariant measure, if it exists, is a Gaussian measure in an appro- 
priate Hilbert space that can be written down explicitly 11, [lit] Assuming 
that the process 

d t z = Az + Q£ 

is ergodic with Gaussian invariant measure /i with mean and covariance 
operator ^A~ 1 Q 2 , then the slow process x converges to the solution of the 
averaged equation 

x = a(x), a(x) = J a(x,y) fi(dy), (34) 

where the integration is over an appropriate Hilbert space. 

When F(-) in ([T]) is given in terms of a symmetric bilinear map, i.e., 
F(v) = B(v,v) the calculation of the vector field that appears in the aver- 
aged equation reduces to the calculation of Gaussian integrals and can be 
performed explicitly. In this case we have 

P c B(x, y) := a(x, y) = D(x, x) + C(x, y) + E(y, y), 



4 The analysis presented in also applies to the case where the fast process is given 
by a semilinear parabolic SPDE. In this more general case, however, it is not possible to 
obtain an explicit formula for the invariant measure of the fast process. 
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where 

N 

D m (x,x) = B k £ m XkXj, 
k,e=i 

N oo 

k=i e=N+i 

oo 

E m (x, y) = ^ BkirnVkyt, m = 1, . . . N, 

k,£=N+l 

and where we used the notation B^ m := (B(ek,ee),e m ) and N := dim (A/") 
denotes the dimension of the null space of A. Then, the fast-slow system (111) 
becomes 

x = D(x,x) + C(x,y) + E(y,y), (35a) 

y = -Ay + b(x,y) + ^ F Q£, (35b) 
e v e 

and the averaged equation for (1351) reads 

x = D(x,x) + E, (36) 

where 

+00 2 

J2 W^ B kkm, m = l,...N. 

k=N+l K 

In the case when the null space is one-dimensional, N — 1, the averaged 
equation becomes 

^ = DX 2 + E, (37) 
at 

2 

with D = Bin an d i? m = J2t=N+i $~Bkki- This equation can be solved in 
closed form: 

x(t) = \ — tan ( VEDt + arctan ( . X ° | ) . 

V £ V Wed) J 

We remark that solutions to (l3Tj) . depending on the choice of the initial 
conditions, do not necessarily exist for all times. We also remark that it is 
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straightforward to consider the case where there is an additional higher order 
linear term (in e) in the equation, i.e. F(v) = B(v, v) + evv. In this case the 
averaged equation fl36l) becomes 

x = D(x, x) + vx + E, 

where x e M N . 

3.3. The Diffusive Time Scale 

We consider the system ([3]) obtained after a diffusive time rescaling to ([T]). 
In order to describe the homogenized equation, we further assume that F(v) 
in (pQ) is of the form 

F(v) = B(v,v) + e 2 uv, (38) 

where £>(•, •) is a symmetric bilinear map satisfying P c B(ek, eu) = 0. 

We recall that the noise does not act directly on the slow variables, 
(Qek,ek) = 0, k = 1 . . . N, where N is the dimension of the null space of 
A. Under appropriate assumptions on the quadratic nonlinearity and on the 
covariance operator of the noise, together with the assumptions on A and Q 
stated earlier in this section, it is possible to prove [7] that the projection of 
the solution to (J3j) onto the null space of A, x :— P c u, converges weakly to 
the solution of the homogenized SDE (the amplitude equation) 

dX = a(X)dt + a(X)dW(t) 1 X(0) = X . (39) 

where the noise is interpreted in the Ito sense and the drift a(x) given by 

a(x) = A^x — B^x, x, x) + vx (40) 

where the linear map A^ : M — >■ jV and the trilinear map B^ : A/" 3 — > M 
are defined by 

A^x = 2B C ((I® S A)-\B S ® S I) + (I®A- 1 B S ) 

+2(5 c ®^" 1 )))(x®Q), (41a) 
B^ = -2B c (x,cA- l B s {x,x)). (41b) 

5 This is essentially the centering condition from homogenization theory, see Equa- 
tion ((291) below. 
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In the above we used the notation B s := P S B and B c := P C B, whereas 
stands for the symmetric tensor product^ and where we have defined 



k=N+l 



k 



The quadratic form associated with the diffusion matrix cr 2 is given by 

+°° „2„2 



(y,<r\x)y)=A £ ll(y,B c (e k ,x)y+ £ ^-— (y,B c (e k ,e^. 

k=N+l k,t=N+l ^ k l > 

(42) 

Furthermore, the fast process can be approximated by an infinite dimensional 
Ornstein-Uhlenbeck process. The precise statement and proof of the above 
results can be found in Q. 

Remark 3.1. The assumption that the 0(e 2 ) term in ( 138|) is linear is needed 
in order to go from (|T1) to ([5]) after rescaling or, equivalently, to (JSj) . If our 
starting point is the already rescaled SPDE then we can apply the results 
from fjj to nonlinearities of the form F(v) = B(v,v) + e 2 h(v) where h(-) 
is an arbitrary nonlinearity . In this case the drift term in the amplitude 
equation (l4"0l) becomes 



ax 



AvoX - B^x.x.x) + J P c h(x,y) fi(dy), (43) 



where n{dy) denotes the invariant measure of the fast OU process. 

When the null space of A is one dimensional and, consequently, the ho- 
mogenized SDE is a scalar equation, it is possible to obtain sharp error 
estimates and to prove convergence in the strong topology. In this case 
Equation fl39|) becomes 

dX = a(X) + a(X)dW, X(0) = (u , e x ) , (44) 



6 Given a Hilbert space H. we denote by r H.® s 'H. its symmetric tensor product. Similarly, 
we use the notation V\ ® s V2 = \{v\ ® V2 + t>2 (8 fi) for the symmetric tensor product of 
two elements and (A (g) s B){x <E> y) = \{Ax ® By + By ® Ax) for the symmetric tensor 
product of two linear operators. Furthermore, we extend the bilinear form B to the tensor 
product space by B(u ® v) = B(u, v). More details can be found in @ Sec. 4]. 
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where 

a{X) = A^X - B^X 3 , a(X) = ^C^ + D^X*. (45) 

In the one dimensional case the formulas for the coefficients that appear in 
the homogenized equation have a simpler form than in the multidimensional 
case. In particular, we have, with B k e m = (B(e k ,ee),e m ): 



oc 



00 ~ ^ + A^Ta7aT' (46a) 



k=2 k,£=2 kl=2 

oo 



^B k \\B\\ k iAax,\ 
Boo = - \ > ( 46b ) 



k=2 K 



oo 



\ ^ zn kml q k q m _ 4tf kn q k (ac„\ 

It is worth mentioning that if we are using a non-orthonormal basis, i.e. a 
basis e k = c k e k , then the coefficients that appear on the right hand side of 
the above equation transform according to 

Bum = -^-B k £ m . (47) 

We also have q k = c k q k . 

Remark 3.2. The formulas for the coefficients that appear in the amplitude 
equation f l5U|) can be also obtained by writing the SPDE §5§ in Fourier space, 
truncating and then using singular perturbation theory-type of techniques for 



the corresponding backward Kolmogorov equation \2u , \21i l. More details on 



this approach can be found in 124} ■ We also remark that, in general, both 
additive as well as multiplicative noise will appear in the amplitude equation, 
although only (degenerate) additive noise is present on the SPDE (CQ). 



4. Numerical Experiments 

In this section we apply our numerical method to several SPDEs and 
report its convergence and performance. We consider here several examples 
of SPDEs with quadratic nonlinearities and check that the theory developed 
in [3] and summarized in Section [3] applies. For all of these examples we can 
derive rigorously the homogenized equation, with explicit formulas for the 
coefficients and therefore, we can present a rigorous numerical study for our 
algorithm and test the effectiveness of the proposed numerical algorithm. 
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4-1- Theoretical considerations 

We will consider variants of the Burgers and the Kuramoto-Shivashinsky 
(KS) equations (with a linear instability term added) in one dimension with 
either Dirichlet or periodic boundary conditions. In particular, we will con- 
sider the singularly perturbed SPDEs (i.e. we have already rescaled to the 
diffusive time scale) 



111 

d t u = —(dl + l)u+-ud x u + vu+-Q£, (48) 
e z e e 

and 



d t u = -(-^ - dt)u + -ud x u + isu + -Q£, (49) 
e z e e 

respectively, where the noise £ is as in Section [3j The operator Q, the co- 
variance operator of the noise, has eigenvalues {qk}kLi and eigenfunctions 
{e^kLx, which are also the eigenfunctions of the differential operator that 
appears in either (148]) or ( 1491) . i.e. the two operators commute. We will con- 
sider these two equations either on [0, it] with Dirichlet boundary conditions 
or on [— 7T, it] with periodic boundary conditions. 

Remark 4.1. For the Burgers nonlinearity and for the boundary conditions 
that we consider it is straightforward to check that the centering condition 
P c B(ek, ejfc) = is satisfied. A more natural equation to consider than ( 1491 ) 
would be the KS equation in the small viscosity regime, i.e. 

d t u = \{-d 2 x - ndi)u + -ud x u + ~Q£, 
e z e e 



where fj, — 1 — v 6 (0, 1). This equation can be rewritten in the form 

d t u = - dt)u + -ud x u + udtu + -Q£. (50) 

e z e e 

The theory presented in fj] and the numerical scheme developed in this paper 
apply to this equation. The application of the numerical method developed in 
this paper to Equation (|50|) and to related models will be presented elsewhere. 
Some recent analytical and numerical results on the behaviour of solutions 



to O50p have been reported in ]3w 
We will use the notation 



A B = {d 2 x + l) and A KS = -d 2 x - d 4 x . 
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It is possible to check that for the above equations and for the chose boundary 
conditions the theory developed in Q and summarized in Section [3] applies. 
Consider first equations ()48p and (|49p on [0,7r] with Dirichlet boundary con- 
ditions. In this case the null space of Ab and Aks is one dimensional: 

Af(A*) = spanj sin(-)}. 

with * being either B or KS. The normalized eigenfunctions of Ab and Aks 
are = y^sin(7rA;). The corresponding eigenvalues are 

Af = k 2 - 1 and Af 5 = & 4 - fc 2 , for fc = 1, 2, . . .. 

Since the null space is one- dimensional, the homogenized equation is a scalar 
SDE. For the nonlinearity B[u, v] = ^d x (uv) it is straightforward to calculate 

B k £m = (B(e k ,e e ),e m ). We have 

B k lm= 1— (|fc + £\Sk+e,m. - \k - £\S\ k -£\,m) , (51) 



2V27T 



oo 

fc=l 



where She denotes the Kronecker delta. We can then use formulas (|46|) to cal- 
culate the formulas that appear in the homogenized equation. Let {— A^} 
of either Ab or Aks with Dirichlet boundary conditions. The homogenized 
equation is given by (144j) that we recall here for convenience 

dX = a(X) + a(X)dW, (52) 

where 

a(X)=A OQ X-^-X 3 , <j(X) = J2 (J^Xt + C^y (53) 
The coefficients that appear in ( ]53"|) can be computed as 



C' 



\ 8 £^ (Afe+i + Afc)AfcAfc +1 

-I 2 2 \ 

16^ A fc A fc+1 (A fc + A fc+1 ) I ' 



7 We use the non- normalized basis e k = sin(7rfc) and use formula ((47j) . 
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In the case where only the second mode is forced with noise, q 2 = <r, qu = 
0, M = 3, . . . then the coefficients become 

In this case only multiplicative noise appears in the homogenized equation 
and it can lead to intermittent behavior of solutions as well as noise induced 



transitions 30 



We will also consider either the Burgers or the KS equation on [—it, it] 
with periodic boundary conditions. In this case the null space of both Ab 
and Aks is two-dimensional and is spanned by 

Af(A*) = span{sin(-), cos(-)}, 

with * being either B or KS. The homogenized equation is given by (I39|) . 
where X = (X±, X 2 ). It consists of a system of two coupled SDEs. We can 
use formulas (j4"U|) and fj42l) . together with the formula for the nonlinearity 
B[u, v] = \d x {uv) to calculate the coefficients that appear in the homoge- 
nized equation. 

4-2. Numerical Experiments 

We shall now apply our numerical algorithm to the model problems , 
(fl9~|) described in Section 14.11 As the behavior of our algorithm is similar for 
the Burgers and the Kuramoto-Shivashinsky equation we will do a thorough 
numerical study on the Burgers equation and comment on the results for the 
Kuramoto-Shivashinsky equation. 

Burgers Equation. We consider equation (1481) on [0, tt] with homogeneous 
Dirichlet boundary conditions. We know from Section 14.11 that, for e suffi- 
ciently small, we have that 

u(-,t) « X(t)sin(-), (55) 

where X(t) is the solution of fl52|) . The function a(X), a(X) in fj4"5l) de- 
pends on Aoo, Coo which for the Burgers equation can be computed using 
formulas (JMD with \ k = k 2 - 1. They read = 0.0026744369, = 
0.00026592835. 

Following the algorithm described in Section ([2]), we look for a solution 
to (SHI) of the form 



M 

u{; t) ~ x(t) sin(-) + Vk(t) sin(/c-), (56) 

k=l 
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substitute the expansion in (T4"8"j) to obtain the a fast-slow system of SDEs as 
described in (|T7|) . Following the algorithm of Section [5] we compute numeri- 
cally the slow variable X n as 

X n+1 = X n + Atd n M + a n M AW n , (57) 

where a M , a M are given by (I23ap and (123b|) . respectively We also consider 
the truncated homogenized problem, i.e., 

dX = a M (X)dt + a M (X)dW, (58) 

where 

d M (X) = A M X-^X 3 , a(X) = ^2 (J-X^ + Cm)^, (59) 

and where A M , Cm, are obtained from (I54al) . (l54bl) with the sums truncated 
at M. 

For numerical comparison we also compute 

■Xn+i,taf = X nM + Ata(X nM ) + a{X nM )AW n , (60) 

-^n+l.hom = X nthom + Ata M (X n)hom ) + aM(X n ^ hom )AW n , (61) 

the Euler-Maruyama approximation of the SDEs fl52l) and fl58l) . respectively. 
The same Brownian path will be used in f|52|) . fl57|) and f )58|) . We emphasize 
that the numerical solutions for fl52|) and fl58|) rely on analytically computed 
homogenized coefficients, whereas for (1571) we implement the multiscale algo- 
rithm of Section El where the coefficients a M ,a M are computed "on the fly" 
and rely on the microsolver (I25ap and (I25b|) . Hence no a-priori analytical 
knowledge of the amplitude equation is required. 

We choose the values of the various parameters entering in the averaging 



process for the computation of a M , a M as suggested in [14J, i.e., K = 1 
5t/e 2 = 0(2~ p ),n T = 0(1), L = 0(2 3p ),L' = 0(p ■ 2 P ). According to [H 
this guarantees (for the case of non-degenerate fast processes) that the error 
is bounded by 2~ p . In our case with a degenerate fast process an error bound 
is still to be established. Here we monitor such convergence numerically. 
More precisely, we set ht = 16, L = 2 3p , L' = p ■ 2 P and monitor the error 
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using 



N 

71=1 

N 



1 N 

TF y^d^Af ~ OM^n.hom)! + I^M ~ (^n,hom) |) (62a) 



^7 E - a(X n , inf )| + \a n M - a(X n>hom )\) , (62b) 



n=l 



for various values of p, where At = T/N and T represent the final time. 
Notice that (I62al) captures the error between (}58j) -the homogenized solution 
of the truncated system-and the numerical solution of the truncated system, 
while f!62bl) . where the index I stands for limit, captures the error between the 
homogenized solution of the limit problem fl52l) and the numerical solution 
of the truncated system. 



2-mode truncation. We set M = 2 in (1561) and substitute the expansion 
in ( 1481) to obtain the following system of equations 

x = vx - -^{xyi + 2/12/2), (63a) 



2e 



(63b) 



( 3 \ 1 / 1 2 x qi . . 

2/i = U-^l2/i--(x2/ 2 --x j + -ei(t), 

2/3 = ( 1/- ^ ) 2/2 + |"(z2/i) + -&(*)■ (63c) 



e 2 y 2e v ai/ e 
The auxiliary process can be derived as explained in Section [2] and reads 

yl = -4ft 1 + -&(*), ( 64a ) 

2/2 = -^2 + -6(t), (64b) 
e z e 

Vl = -^yl-k(xyl-%r), (64c) 



e 2^ e 2 y - 2 , 
8 3 

2/2 = -#vl + 2# x yi' ( 64d ) 

We apply the algorithm of Section [2] to get a numerical approximation of the 
homogenised problem corresponding to fl63l) . The final time is T = 1 and 
N = 10, which corresponds to macro time-step of size At = 0.1. The macro 
solver for the method is given by 057p . As mentioned above, we compare our 
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results with (1601) and (IBTj) . The unknown coefficients A3, C3 in (l5T)j) can be 
computed using f )54ap and ( 154bj) . where the sum is truncated at M + 1 = 3 
and read A 3 = 0.003735726834, C 3 = 0.0002593873518. 



10° 
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Figure 1: Numerical convergence for 2-mode truncation. On the horizontal axis we monitor 
the accuracy of the micro-time step and on the vertical axis we measure the error as given 
by (|62a]) and ((62b)) with M = 2. 

We observe in Figure [1] that we get numerically the expected order of 
convergence corresponding to 5t/e 2 = (D(2~ p ). Furthermore, as the micro 
time-step becomes smaller, the numerical scheme gets closer to ( IBTl) and 
slightly deviates from ( 1521) . This is expected as the numerical solution is 
not converging to that latter solution. We observe nevertheless that with 
only two fast modes, the numerical scheme already captures quite well the 
effective behavior of the slow variable of the infinite dimensional system. 

We also illustrate the time evolution of one trajectory comparing over the 
time < t < T with T = 10, the Euler-Maruyama method for the amplitude 
equation f[6"Uj) . the homogenized equation ( IBTj) and the macro solver (I5T|) . 
The same Brownian path is used for generating the three trajectories and as 
well as the same macro time step. We perform this comparison for increasing 
accuracy of the micro solver used to recover the macro data, namely, 5t/e 2 = 
0{2~ p ), p = 3,4, 5. We see in Figure [2] that the trajectory for the amplitude 
equation and the homogenized equation coincide, while the macro solver gets 
closer to the true dynamics as we refine the micro time step. For the same 
trajectory we also give a space-time plot for the approximation of the original 
SPDE u(-,t) w X{t) sin(-), with X(t) solution of the amplitude equation, the 
homogenized equation or the macro solver. Again we see that the numerical 
method captures the right behavior of the solution. 
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Figure 2: Euler-Maruyama methods (|57j) (solution denoted X n ), (jFTj) (solution denoted 
l„ji 0m ) and ([SO]) (solution denoted X n ^ n f) for three paths (left p = 3 for X„, middle 
p = 4 for X„, right p = 5 for We use 2-mode truncation for (f5Tf and ([STjl . 



3.5 3 " 5 \ 3 '•' 




o o oo o o 



Figure 3: Approximation (|55|) of the solution u(x,t) of the SPDE; u(-,t) ~ X n (t) sin(7r-) 
(left figure p = 3), — -^„,hom(0 sin(7r-) (middle figure) and u(-,t) ~ X„ j n f(i) sin(7r-) 

(right figure). 



3-mode truncation. We set M = 3 in ( l56|) and obtain the following system 
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of equations 



x = ux - — (xyi + 2/12/2 + 2/22/3) , (65a) 




The auxiliary process can be computed similarly as for the 3-mode trunca- 
tion. We perform the same set of numerical experiments as for the 3-mode 
truncation, reported in Figure 14.21 Similar behavior than previously noted 
can be observed. Observe that the discrepancy between the numerical scheme 
and fl52|) gets smaller. This is expected as with additional modes, the ho- 
mogenized equation f )58|) (that we aim at approximating with our multiscale 
scheme) gets closer to ( |52l . 
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Figure 4: Numerical convergence for 3-mode truncation. On the horizontal axis we monitor 
the accuracy of the micro-time step and on the vertical axis we measure the error as given 
by (|62a) and ([62b]) with M = 3. 

4-mode truncation. We set M = 4 in ( 1561) and apply the similar procedure 
as previously. For the sake of brevity, we do not write the system of equations 
in this case and just report the numerical convergence. 

We see in Figure I4.2I a similar behavior of our numerical scheme as ob- 
served previously. We again notice that the discrepancy between the numer- 
ical scheme and ( |52i) is smaller than for lower order truncation. Notice that 
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the first numerical result reported is obtained for St/e 2 = 2~ p with p = 4. 
This is due to stability issues with the Euler-Maruyama scheme for the fast 
process. As the linear term in the equation for the fifth mode y$ is (y — ) 7/5, 
the stability restriction 24St/e 2 < 2 implies St/e 2 < 1/12 and thus the time 
step 5t/e 2 = 1/8 is unstable. 
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Figure 5: Numerical convergence for 4-mode truncation. On the horizontal axis we monitor 
the accuracy of the micro-time step and on the vertical axis we measure the error as given 
by (|62a| and ((62b)) with M = 4. 




The Kuramoto-Shivashinsky equation. . The equations for the M-mode trun- 
cation of the Kuramoto-Shivashinsky equation are very similar to the ones 
for the Burgers equation and will not be presented here. The only difference 
is that the fast process is more dissipative than for the Burgers equation, due 
to the stronger dissipativity of the operator Aks compared to As- As the 
results of the numerical experiments for the KS equation are very similar to 
the results reported in this section for the Burgers equation, they will not be 
presented here. 

5. Conclusions and Further Work 

We have presented a new numerical method for the efficient and accurate 
solution of stochastic partial differential equations with multiple scales. The 
new numerical scheme is based on a combination of a spectral method with 
the HMM methodology and has been tested on SPDEs with quadratic non- 
linearities for which a rigorous homogenization theory exists. This enables 
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us to check the performance of our method. The numerical experiments pre- 
sented in this paper suggest that the new method performs well and allows 
to solve accurately multiscale SPDEs by solving a low dimensional fast-slow 
system of SDEs. The method is suitable for infinite dimensional stochastic 
systems for which there is clear separation of scales, and for which a low 
dimensional homogenised (or averaged) equation for the slow modes exists. 

There are still many questions that are left open. First, the rigorous 
analysis of the proposed method and a careful study of the convergence 
and stability properties of the proposed method remains to be done. In 
addition, the optimisation of the proposed method by tuning appropriately 
the parameters of the algorithm has not been performed yet. This appears 
to be an open problem even when the HMM methodology is applied to finite 



dimensional fast / slow systems of SDEs [22 . 

The proposed numerical algorithm could be used to study in detail the 
qualitative and quantitative properties of solutions to SPDEs with quadratic 
nonlinearities, since SPDEs of this form exhibit very rich dynamical be- 
haviour. Furthermore, we would like to apply the numerical algorithm to 
more general classes (and systems) of semilinear SPDEs, for which an aver- 
aged or homogenised equation is known to exist. Examples include systems 



of reaction/diffusion equations that were considered in [10] as well as the 
Swift-Hohenberg SPDE @. 

In our algorithm, we did not make use of the fact that the form of the 
amplitude equation (i.e. a Landau equation with additive and multiplicative 
noise) is known. Knowledge of the functional form of the coefficients that 
appear in the homogenised or averaged equation can be used in order to 
simplify the numerical algorithm. The stochastic Landau equation appears as 
the amplitude equation for several infinite dimensional stochastic dynamical 
systems, not only for SPDEs with quadratic nonlinearities, e.g. [ij]. Thus, the 
algorithm proposed in this paper, could be modified to develop an efficient 
method for studying the dynamics of SPDEs near bifurcation. All these 
topics are currently under investigation. 
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